Relevance of complete Coulomb interaction matrix for the Kondo problem: 

Co impurity on Cu(lll) 
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The electronic structure of a prototype Kondo system, a cobalt impurity in a copper host is 
calculated with accurate taking into account of correlation effects on the Co atom. Using the 
recently developed continuous-time QMC technique, it is possible to describe the Kondo resonance 
with a complete four-index Coulomb interaction matrix. This opens a way for completely first- 
principle calculations of the Kondo temperature. We have demonstrated that a standard practice 
of using a truncated Hubbard Hamiltonian to consider the Kondo physics can be quantitatively 
inadequate. 
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Introduction 

Scanning tunneling microscopy (STM) has become one 
of the most basic tools for the manipulation of mat- 
ter at the atomic scale. Although this experimental 
technique has reached maturity, the detailed theoretical 
understanding of experimental data is still incomplete 
and/or contradictory. One of the most famous examples 
of atomic manipulation is associated with the surface 
Kondo effect observed when transition metal ions (like 
Co) are placed on a metallic surface (such as Cu (lll))^ 2 -. 
The surface Kondo effect is the basis for the observation 
of surprising phenomena like quantum mirages^, and has 
attracted a lot of attention and interest in the last few 
years. Early interpretations of these observations were 
based on the assumption that only surface states of Cu 
(111) are involved in the scattering of electron waves by 
the Co adatoms^ 5 ^. However, later experiments with Co 
atoms on the Cu (100) surface (that does not have any 
surface state)£, or in Cu (111) but close to atomic sur- 
face steps (that affect the surface states )£ have indicated 
that bulk rather than surface states are responsible for 
the Kondo effect in these situations. The latter can be 
important for fine tuning of surface electronic structure, 
with potential applications to nanotcchnology. A recent 
study of CoCu„ clusters on Cu (111) demonstrated this 
tunablilty by atomic manipulation and showed that each 
atom in the vicinity of the magnetic impurity matters for 
determining the Kondo effect^. Moreover, the relevance 
of the Kondo effect for the electronic structure of metal 
surfaces themselves was demonstrated by the discovery 
of a sharp density of states peak on the Cr (001) sur- 
face and its possible interpretation as an orbital Kondo 
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At the same time, when calculating the Kondo tem- 
peratures for real electronic structures a mapping onto 
one-orbital Anderson impurity modeU^ was used. The re- 
alistic atomic geometry of Kondo systems plays a crucial 



role in complex electronic properies 9 -^ and it is, a pri- 
ori, not obvious that a one-orbital Anderson impurity ap- 
proach is sufficient: even the two-orbital Anderson model 
demonstrates Kondo physics essentially different from 
the single-orbital onei£. A recent theoretical investiga- 
tion of Fe impurities in gold and silver showed that the 
proper Kondo model corresponds to a S=3/2 spin stated. 
A realistic, multi-band consideration of correlation ef- 
fects in specific solids is possible in the framework of the 
Local Density Approximation + Dynamical Mean-Field 
Theory (LDA+DMFT) approach (for review, see Ref. 
[It]). However, formally accurate Quantum Monte Carlo 
(QMC) calculations are always done with taking into ac- 
count only the diagonal part of Coulomb interactio n 18 : 19 , 
even with realistic hybridization functions obtained in 
the LDA. This approximation is, strictly speaking, un- 
controllable. At the same time, approximate schemes 
working with the complete Coulomb interaction matrix, 
such as the perturbative scheme 2 ^ which is frequently 
used to calculate electronic structure of transition met- 
als and alloy a 21 ' 22 are not sufficient to reproduce so subtle 
correlation features like the Kondo effect, properly. As 
for the exact diagonalizatio n ' or numerical renormal- 
ization group 11 ' 15 ' 23 methods they are hardly applicable, 
due to computational problems, for more than two or- 
bitals per impurity. 

The recent progress in continuous time QMC scheme 
(CT-QMC) 24 ' 25 makes it perspective to treat the com- 
plicated Kondo systems 2 ^. Here we will apply this 
method to calculate Kondo temperatures as well as spec- 
tral functions for the case of a Co impurity in bulk Cu, 
in a Cu (111) surface and on top of a Cu (111) sur- 
face. In contrast with all previous calculations we will 
work with an accurate complete Coulomb interaction U - 
matrix for correlated d orbitals. The latter can be cal- 
culated from first principles in a parameter-free way by 
the GW technique 2 ^ so this approach is completely ab 
initio. Moreover, the CT-QMC method allows to work, 
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without any essential difficulties, even with the rigor- 
ous frequency-dependent {/-matrix. As the first step, we 
present calculations for the static [/-matrix, but this re- 
striction is purely technical and can be relatively easily 
removed in the future, with a growth of available com- 
puter resources. 



here e n k is the energy spectrum and V'nk is the corre- 
sponding wave function of our system (metal host with 
magnetic impurity), described by di localized orbitals, 
and Uijki is the Coulomb interaction matrix element: 



U m = (hh\V{ 2 e \k 2 h) 



(3) 



I. MULTI-ORBITAL CT-QMC FORMALISM 

The multi-orbital impurity problem with a general U - 
matrix is described by the effective action 
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where i,j,k,l are orbital indices, and a, a' are spin in- 
dices, Qij is the local non-interacting Green function for 
correlated orbitals obtained from the Density Functional 
Theory (DFT) with the help of optimal projection oper- 
ator to the impurity d-states: 
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here i\ = di (ri) is local orthogonal wave function for 
correlated orbitals and is screened spin-independent 
Coulomb interaction between electrons at the coordi- 
nates ri and r 2 . We used standard quasiatomic LDA+U 
paramctrization of Coulomb matrix for d-electron via ef- 
fective Slater integarls or average Coulomb parameter U 
and exchange parameter J as described in RefUg). We 
choose the orbital basis related to spherical harmonics to 
be sure that magnetic orbital quantum numbers in Uijki 
matrix are satisfied the following sum rule: i + j = k + 1. 
In this case we will get rid of so-called three-site terms like 
Uikki with i ^ I which turns out to result in a strong sign 
problem in QMC calculations with real spherical harmon- 
ics. 

Following the general CT-QMC scheme^ we expand 
the partition function around the Gaussian part of our 
multiorbital action Eq.([T]) which gives the fermionic de- 
terminant over the non-interacting Green functions with 
the rank In: 
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In order to minimize the number of different interac- 
tion vertices we group different matrix elements of the 
multiorbital Coulomb interactions which have a similar 
structure of fermionic operators. Since the Uijki matrix 
elements are spin independent, one should look over all 
possible combinations of orbital and spin indices, to gen- 
erate all terms for the interaction in the action Eq.([T]). 
Some combinations can violate the Pauli principle and 
should be removed. For CT-QMC algorithm it is useful 
to represent the interaction Hamiltonian in the following 
form: l}ij k ic\ a ci a c^ al c ka ' . 

The interaction terms can be transformed to the de- 
sired form, depending on relations between spin and or- 
bital indices: 

(i) if a ^ a' . we can just commute c\ a and Cka' and 
then ciu and c^ a , . Another combination of indices, that 
allows the same commutation, is the following: a = a', 
i ^ j and k / I (the latter two are following from the 
Pauli principle), and also j ^ I. These terms we can 
transform to the following desirable representation: 
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(5) 



(ii) in the case when a = a' and j = I we can commute 
Cko' and ct^,, since in this case i ^ j and k ^ I due to 
the Pauli principle: 



Tjint2 _ _tj . J „ J „ 



(6) 



After generating all this terms it is useful to collect 
and symmetrize all the terms with identical and equiv- 
alent (i.e. Ui ]k ic\ a c 3<J c{ a ,ci a , and Vmi^c^^c^) 
quantum numbers. 

In order to reduce the fermionic sign problem we intro- 
duce additional parameters, a, to optimize the splitting 
of the Gaussian and interaction parts of the action Eq. JT]) 
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One can see, that the first item in ([7]) on Matsubara 
frequences corresponds to bare Green's function 

Gij 1 = {iUn + fj) Sij - Aij(u> n ) (8) 

where A is the hybridization matrix. The second term is 
just a constant which we can absorb to the new chemical 
potential fl. Therefore we can rewrite the bare Green 
function in the following matrix form: 

Q' 1 = (iuJ n + p) 1 - A, (9) 

The optimal choice of parameters afj would lead to ef- 
fective reduction of interaction terms in the action Eq. 
(|7|) and therefore minimization of average perturbation 
order in Eq. (j4|). 

Note that relation between Q and Q can be represented 
from Eq. (|7|) in the following spin and orbital matrix 
form: 

Q- 1 = Q- 1 - (afj) . (10) 

Here we used the fact that Uukj = Uujk following from 
the definition of the Coulomb matrix elements ©. 

We also need to minimize the fcrmionic sign problem 
which finally leads us to such expression for diagonal al- 
pha parameters 

a« + a&=a, (11) 

corresponding to to the following interaction fields 
Uijjini a rij a i . The a has to be found iteratively in order to 
get a proper occupation number of correlated electrons. 
In the case of half-filled one-band Hubbard model a — 1 
leads to the correct chemical potential shift of the ¥ and 
average a — i which corresponds to the Hartree-Fock 
substraction. For non-diagonal alpha's which correspond 
to the fields of general form Uijkic\ a cic,c^ al Ckcy' , where 
i =/= I and j k we find the following optimal condition: 

o% + c$ = (12) 

Since we symmetrize the interaction U matrix it is nec- 
essary to extend the definition of the a matrix to keep 
all the terms in the interaction part of initial action ( the 
last item in Eq. (J7J)). It can be done in the following 
wa y 24 i 34 : for every Uijki field in 50% of updates we de- 
liver the a parameters as 

a = ctdiag, & = a — adiag, and in another 50% as 



a d = a — adiag, = adiag for the case of i = I 
and j = k. For non-diagonal fields, i.e. i ^ I and 
j 7^ k a a = a n d, a jfc = —a n d, with 50% probability 
and a %l — —a n d, a^ k — a n d otherwise. It was found that 
the sign problem is eliminated in the case when adiag < 
and a > 1 for occupancy n > | per state and adi ag > 0, 
a < 1 otherwise. The optimal choice of |a^ a g| param- 
eter is few percent of \a\ to keep minimal average per- 
turbation order. Another problem is a proper choice of 
non-diagonal a n d parameter. It is easy to see that a n d 
is proportional to acceptance probability of non-diagonal 
field in the case where corresponding bare Green func- 
tion Qjk = 0. Since these processes are unphysical, the 
natural choice is a n d = 0. But it leads to division by zero 
in the updating the inverse Green function matrij*21. On 
the other hand increasing the a n d parameter causes a 
sign problem. We find a reasonable choice of a n d to be 
on the order of 10 -4 . Moreover for some special cases 
like the atomic limit, where Qmm(i~) is constant, a small 
noise should be added to all the a parameters to avoid 
numerical divergency. 



II. RESULTS 

The Co-Cu system is treated as five-orbital impurity 
model representing 3d electronic shell of the cobalt atom 
hybridized with a bath of a conduction Cu electrons. 
The bath Green function was obtained using the first- 
principle density-functional theory within the supercell 
approach. For Co impurity atoms in the bulk as well as in 
/ on the Cu (111) surface the bath Green functions were 
obtained using the Vienna- Ab-initio simulation package 
(VASP) 35 ' 36 using the projector augmented wave (PAW) 
basis sets^i. The density functional calculation for cobalt 
impurity in the bulk was carried out using CoCu63 super- 
cell structure with lattice constant corresponding to pure 
copper. The surfaces were modelled by supercells of Cu 
(111) slabs containing 5 Cu layers with 2x2 and 3x4 
lateral extension for Co in and on the surface, respec- 
tively. The PAW basis naturally provides the projectors 
(di\ipnk) required in Eq. @. In using these PAW pro- 
jectors, directly, we employ here the same representation 
of localized orbitals as used within the LDA+U-schemc 
implemented in the VASP-code itself or as discussed in 
the context of LDA+DMFT in Ref. [M 

For the problem of a single Co impurity in a bulk cop- 
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FIG. 1: (color online) Comparison with ED in the atomic 
limit (without hybridisation to the bath of free electrons). 
Main graph: U = 1 eV, J = 0.4 eV /3 = 2 eV" 1 , for 5-orbital 
impurity at half-filling; inset: (7 = 2 eV, J = 0.7 eV, j3 = 3.7 
eV -1 , for 5-orbital impurity with 8 electrons. 



per matrix the basis set of spherical harmonics Yi m is 
used. In this basis the interaction part of the hamil- 
tonian contains only terms of the following form: diag- 
onal density-density like H^ 9 = Uijjin ia nj a i , where 

n ia = c\ a Cia and non-diagonal = Uijkl^Cia'^c-ic, 
where i ^ j and k ^ I. The Coulomb matrix for the d- 
electron shell in the basis of complex harmonics contains 
45 non-equivalent diagonal terms. Non-diagonal terms 
can be further classified to a spin- flips, where i = I, 
j = k and the most general four-orbitals interactions, 
where this condition is not fulfilled. Notice, that pair- 
hopping terms (i = k, j = I) are restricted by symmetry 
in this basis. In description of e?-electron shell we have to 
involve 20 non-equivalent spin-flips and 64 terms of the 
most general form. 

To find the effects, caused by non-diagonal terms, 
we used two different interaction Hamiltonian. First, 
interaction with only diagonal terms was used. In this 
case there is no sign problem. Then, the complete 
Coulomb interaction matrix of the 3d-electron shell of 
the cobalt atom with 129 terms was included. 

As a benchmark we use impurity problem in the atomic 
limit, since it can be compared with the result of ex- 
act diagonalization (ED) method. The results imaginary 
time Green function for 5-orbitals model with different 
chemical potentials corresponding to the d 5 and d 8 con- 
figurations are shown in the Fig. [1] in comparison with 
ED results. The significant difference between density- 
density (diagonal) interaction and the full vertex can be 
found both at half-field case with relatively high tem- 
perature with the [7 = 1 eV, J = 0.4 eV, /3 = 2 eV' 1 
and at non-symmetric case even for lower temperature. 
Note that in the d 8 and d 7 cases the many-body ground- 



FIG. 2: (color online) Histograms of Monte-Carlo distribu- 
tions for average perturbation order. Main graph: (7 = 4 eV, 
J = 0.7 eV, j3 = 10 eV -1 for 5-orbital impurity coupled to 
realistic Cu-bath with 7 electrons; inset: (7 = 4 eV, J = 0.7 
eV, beta = 1 eV -1 in the case of 5-orbital impurity model, 
coupled to semi-elliptical bath with bandwidth W = 0.5 eV 
at the half-filling 



state have different symmetry for diagonal interactions 
and non-diagonal full vertex. The results for d 8 configu- 
ration with the interaction parameters U = 2 eV, J = 0.7 
eV, (3 = 3.7 eV^ 1 are shown in the insert to the Fig. Q] 
The difference between Green function of the interacting 
system with full Coulomb interaction and density-density 
one is visible on the G(t). We find a very good agreement 
between CT-QMC results and ED solution. 

In the inset of Fig. [2] we show the distribution of non- 
diagonal terms, i.e. the contribution of Coulomb fields of 
the form (JSJ) to the resulting Green function. The zero 
entry of this histogram counts the number of steps when 
all the fields contributing to the fermionic determinant 
(|4"|) were of density-density type. The entry with index 2 
show us the number of steps where the average ([J) was 
containing two spin-flip type fields J5]). Such situation 
takes place, for example, when one Coulomb field rep- 
resenting spin-flip cj|Cj^cj|Cjj process was used to con- 
struct the determinant. 

One can see in the insert of the Fig. [21 that only 
even orders of interaction histogram have large accep- 
tance probability at high temperature and even the tenth 
order in non-diagonal interactions has non-zero contribu- 
tion. The 3-rd and 5-th order contributions exists due to 
the finite a nd parameter. 

Typical distribution of the perturbation order (5- 
orbital AIM with 7 electrons, U = 4 eV , J = 0.7 eV, 
[5 = 10 eV^ 1 ) is shown in Fig. [2 main plot. Dash line de- 
notes the perturbation order during accepted steps that 
involved non-diagonal fields. The coincidence of distribu- 
tions maxima of both histograms demonstrate that the 
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FIG. 3: (color online) Total DOS of 3d orbital of Co atom 
embedded in Cu matrix. Model parameters: U — 4 eV, J 
0.7 eV, P — 10 eV -1 for 5-orbital impurity with 7 electrons. 



acceptance rate mostly depends on diagonal interactions. 

For many-body calculations of the Co impurity in the 
Cu matrix we need to find the effective d-orbital chemi- 
cal potential which defines the number of 3d-electrons of 
cobalt. The particular electronic configuration of a Co 
atom in a copper matrix is unknown, but the DFT re- 
sults (rid — 7.3) give us an evidence that it is close to 
d 7 configuration. Therefore we performed all impurity 
calculations for cobalt d 7 configuration. 

The results of the CT-QMC calculations for U = 4 
eV and J = 0.7 eV are presented in Fig|3] compared to 
the bare impurity density of states for cobalt impurity 
in the bulk. There is a pronounced difference between 
Kondo-like resonance near the Fermi level. In the case of 
full U- vertex it becomes more narrow and located much 
closer to the Fermi level. The sign problem for realis- 
tic five-band model depends crucially on the symmetry 
of coulomb Vertex Uijki and magnitude of non-diagonal 
terms in the bath Green functions Gij . The most serious 
problem is related with non-diagonal terms of U-matrix, 
therefore we use a basis of complex spherical harmonics. 
In this case there is no so-called three-cite terms or cor- 
related hopping, e.g. Uikki- On the other hand, in this 
basis, the bath Green-function matrix Qij for d-electrons 
has two non-diagonal elements in the bulk of cubic crys- 
tals and much more on the surface and in the first layer. 
Moreover there arc lot of small four-site terms Uijki which 
result in a large sign problem for surface-adatom calcu- 
lations. The sign problem for a Co impurity in the bulk 
is not large and average sign is between 0.90 and 0.97 
depends on the simulation temperature. 

In the case of non-diagonal interaction we used so- 
called cluster steps which correspond to complex Monte- 
Carlo updates with more than one additional interaction 
field. This scheme became essential for spin-flip like inter- 
action or more general U-vertex which can contribute to 



FIG. 4: (color online)Total DOS of 3d orbital of Co atom 
embedded in the bulk of Cu, into 1-st layer and Co-adatom on 
the Cu(lll) surface. Model parameters: U = 4 eV, J = 0.7 
eV, p = 10 eV" 1 for 5-orbital impurity with 7 electrons. 



the Green- function only in the second or higher order "di- 
agammatic" expansion and this can let the Monte-Carlo 
process to explore all the phase space. We note that prob- 
ability of non-diagonal terms drastically decrease with 
increasing the hybridization to the bath. Nevertheless, 
at least for three-band benchmarks we found remarkable 
effect of the spin-flip terms if the bath Green function has 
peaks in the vicinity of the Fermi level on the distance of 
the order of J. 



We estimated the rcnormalization factor Z = (1 — 
dZ/dE)' 1 for U = 4.5 eV, J = 0.7 eV and = 10 eV" 1 
and find Z t „ = 0.5 and Z e = 0.4 which shows the rea- 
sonable strong interaction of Co d-electrons. We estimate 
the Kondo temperature (Tk) using the temperature de- 
pendence of FWHM for resonance near Fermi level. Since 
our simulation temperature is very high compared to Tk 
we can get only order of magnitude of Tk = 0.1 eV, 
which is reasonable for Co-impurity systems. 

We also performed the CT-QMC calculation of cobalt 
impurity on the surface of Cu(lll) and embedded into 
the first copper layer. In contrast to the bulk system 
the surface one has a large sign problem, related with 
the relativelly large non-diagonal elements of the bath 
Green functions. Although changing of the sign is a very 
rare event (less then 0.03% of the accepted steps) and we 
used a simple constrained sign calculations. Comparison 
of the different spectral functions for the bulk, surface 
and fist-layer cobalt impurity is presented in Fig. 3J One 
can see clearly the change of the Kondo resonance width 
as a function of reduced dimensionality. 
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CONCLUSIONS 

In conclusion, we perform the continuous time QMC 
calculation of realistic 5-orbital Co impurity in cop- 
per and discuss the relevance of non-diagonal part of 
Coulomb matrix in the Kondo problem. Comparing Figs. 
[3] and 5] we find that non-density-density terms in the 
Coulomb vertex are required to obtain quantitative pre- 
dictions of spectral functions and related properties. The 
position of the Hubbard peaks and the Kondo peak is 
markedly changed by spin-flips and other non-diagonal 
terms of the Coulomb vertex. Thus, obtaining sensitive 
observables like Kondo temperatures quantitatively re- 
quires accounting for these terms. On the other hand 
hybridization effects like bringing the Co impurity from 
bulk to the surface and having it on top of the surface can 



be quite drastic: As Fig [4] shows, the sharpening of the 
Kondo resonance and the shifting of the Hubbard bands 
is much stronger when going from bulk to the surface 
then on switching on the non-diagonal part of Coulomb 
matrix. Only the qualitative overall shape the DOS and 
its response to strong hybridization changes are well de- 
scribed by density-density type terms of the Coulomb 
vertex. 
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